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Abstract. The DGMRES method for solving Drazin-inverse solution of sin- 
gular linear systems is generally used with restarting. But the restarting often 
slows down the convergence and DGMRES often stagnates. We show that 
adding some eigenvectors to the subspace can improve the convergence just like 
the method proposed by R.Morgan in [R.Morgan, A restarted GMRES method 
augmented with eigenvectors, SIAM J. Matrix Anal.Appl. 16 (1995)1154-1171]. 
We derive the implementation of this method and present some numerical ex- 
amples to show the advantages of this method. 
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1 Introduction 

DGMRES is a new iterative method for computing the Drazin-inverse solution 
of linear systems [3 . Consider the linear system 

Ax = b (1) 

where A 6 C™ x ™ is a singular matrix. We recall that the Drazin-inverse solution 
of ([1} is the vector A D b, where A D is the Drazin-inverse of the singular matrix 
A. 

In [2], A.Sidi proposed a general approach to Krylov subspace methods for 
computing Drazin-inverse solution. In that paper the authors do not put any 
restriction on A, that is, A is non-hcrmitian or hermitian, index of A is arbi- 
trary and the spectrum of A can have any shape. In [3J, A.Sidi gave one of 
the Krylov subspace method named DGMRES, which is a GMRES-like algo- 
rithm. Like GMRES, in practical use, we often propose restarted DGMRES 
which denoted by DGMRES (m). DGMRES (m) is an economical computing 
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and storgewise method for the Drazin-inverse solution. But restarting slows 
down the convergence and often stagnates (see @],[5]). In the present paper we 
propose a DGMRES method augmented with eigenvectors, which can acceler- 
ate the convergence and overcome the stagnation. Classical GMRES augmented 
with eigenvectors was studied by R.Morgan (see [T]). Now we derive a DGM- 
RES augmented with eigenvectors similarly. We give the convergence analysis 
of DGMRES augmented with eigenvectors which shows this method is more 
effective than DGMRES without augment. Then some numerical experiments 
are presented to show the convergence rate of DGMRES augmented with eigen- 
vectors is remarkably improved, especially when the matrix has a very large or 
small nonzero eigenvalue. 

The paper is organized as follows. In section 2, we will give a brief review 
of DGMRES. In section 3, we obtain the convergence analysis of DGMRES 
augmented with eigenvectors and derive the algorithm. In section 4, we present 
some numerical examples to compare the DGMRES with the DGMRES aug- 
mented with eigenvectors. 



2 DGMRES 

Throughout this paper, we suppose the index of A is known. The index of A, 
denoted by ind(A), is the size of the largest Jordan block corresponding to the 
zero eigenvalue of A. DGMRES is a Krylov subspace methods for computing 
the Drazin-inverse solution A D b. For more details we refer the readers to 

We start with a initial vector xq and the method we interested in is to 
generate a sequence of vectors x\, ■ ■ ■ , which satisfies 

x m = xq + q m -x(A)ro, r —: b - Ax , 

m — a 

where q m -i(X) = £ CiA"^ 1 (a := ind(A)) . Then 

r m :=b- Ax m = p m (A)r (2) 



m — a 



where p m (\) = 1 - Aq m _i(A) = 1 - £ c t A a+l . 

i=l 

The Krylov subspace we will use is 

K m (A; A a r ) = S pan{A a r , A a+1 r , ■ ■ ■ , A"^} 
. The vector x m produced by DGMRES satisfies 

M«r m ||= mm \\A a (b - Ax)\\. (3) 

x£xo+K m (A;A a r ) 

Now we give restarting DGMRES Algorithm DGMRES(m). 
Algorithm 1 DGMRES (m) 

1. Choose initial guess xq, and compute = b — Axq and A a ro; 
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2. Compute (3 = \\A a r \\ and set v\ = /3~ 1 A a r ; 

3. Orthogonalize the vectors: A a rQ, A a+1 r$, ■ ■ ■ , A m ro by Arnoldi- Gram- Schmidt 
process. 

For j = 1, 2, • • ■ do 
For i = 1, 2, • • • do 
Compute hi t j = (AVj,Vi). 

^ j 
Compute Vj — Avj — hijVi. 

Set // , . i . = \\vi\\ and v J+1 — Vj/hj+ij. 

4. For k = 1, 2, • • • „ /orm tfie matrices V k E C nxk and H k G C( fe+1 ) xfe . 





v k = 




■ ■ -,vk\; 




( h n 
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V o 




frfc+l.fc / 



5. For m = a + 1, ■ • /orm £/ie matrix H m = H m H m -i ■ ■ ■ H m — a . 

6. Compute the QR factorization of H m : H m = Q m R m , where R m is upper 
triangular. 

7. Solve the system R m ym = P(Qm e i)> where e\ = (1, 0, • • ■ , 0) T . 

8. Compute x m = x a + V m - a y m . 

9. Restart: if r rn = b — Ax m satisfied with the residual norm then stop, else let 
%q = %m cmd go to 2. 

In [H|3], the convergence analysis of DGMRES is given as follows. 

Lemma 2 Denote the spectrum of A by o~(A) and choose £1 to be a closed 

domain in the A— plane that contains a(A)\{0} but not A = ; such that its 
boundary is twice differentiable with respect to arc-length. Denote by the 
conformal mapping of the exterior of Q onto the exterior of the unit disk {w : 
\w\ > 1}. Then the vector x m extracted from K m (A; A a ro) satisfies 

\\A a r m \\<Km a+2 ^p m , (4) 

for all m, where K is a positive constant independent of m, k = max{kj : kj = 
ind(A-Xj),Xj G a (A) \{0}} and p = |1/$(0)| < 1. 

We know the x' m s produced by DGMRES is 

m — a 

x m = x a + ^2 c i A a+i - 1 r = x +p(A)r . 
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Let Zi, Z2, ■ • ■ , zi is a set of linear independent eigenvectors corresponding to 
eigenvalues Ai, A2, • • • , Aj. We add some new vectors 2j> + 2, ■ • • , z* n such that 

I n 

zi, ■ ■ -zi,zi + i, ■ ■ ■ ,z n is a basis in C". Then r = b- Ax Q = Pi*i + Pi% 

i=l i=l+l 

and 

r m = b- Ax m = b- A(x Q + p(A)r Q ) 

l n 

= r - Ap(A)[J2PiZi + f 3 ^ 
j=i i=i+i 

l n 

= r - ^2\ l p(\ l )z l - £ (3 t Ap(A)z, 

i=l i=l+l 
I n I n 

= ^PiZi+ £ PiZi - ^2\ l p(X l )z l - £ p t Ap(A)z t 

i=l i=l+l i=l i=l+l 

I n 

= £(ftz 2 - Xip(Xi)zi) + Wi ~ PiAp{A))%. 

4 = 1 1 = 1 + 1 



Thus 



I n 

A a r m = ^[ftA^ - A? +1 p(Aj)*i] + [PiK%-PiP(A)A a+1 Zi}. 

i=l i=l+l 

From Lemma 2, we know 

I n 

|| £[ftA^ - X<t +1 P {X t ) Zl ] + £ [PiXIzi ~ p tP {A)A a+1 z^ || (5) 

i=l i=l+l 

< Km a+2 &-Vp m . 



3 DGMRES augmented with eigenvectors 

In order to accelerate the convergence of restarting GMRES, Morgan suggested 
some eigenvectors corresponding to a few of the smallest eigenvalues add to the 
Krylov subspace for GMRES. Then the convergence can be much faster. This 
method can be used in the process of DGMRES. 

Let k be the number of the eigenvectors added to the subspace. Let r be the 
initial residual vector and a = ind(A). K m := span{A a ro, A a+1 r , • • • , A m ~ 1 r } 
and K m . k = span{A a r , A a+1 r Q , ■■■ , 4 m_1 r , Zi,z 2 , • • • , z k }, where zi,z 2 , ■•■ ,z k 
are the eigenvectors added to K m . The approximated Drazin- inverse solution 
will be extracted from K m ^ k . That is 

k m — a 

x m = x + £ a z Zi + £ c l A a+l ^ 1 r 
i=i i=i 
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. which satisfies 



\\A a r m \\ = min \\A a (b - Ax)\\. (6) 

xex +K mik (A;A°-r ) 



Suppose zi, Z2, ■ ■ ■ , Zk, • ■ ■ , zi is the set of linear independent eigenvectors and 
z\ , ■ ■ ■ , Zk, ■ ■ ■ , zi, zi+i, ■ ■ ■ , z n is a basis of C™ . So we can write 



r ° = E ^ lZl + E @ iZi + E @ iZi ' 



i=k+l i=l+l 



and 



= + E aiZi + ^(^) r ° 



k 

i=l 

= x + '^2a i Zi+ p(A) PiZi+ E $ iZ ^ ■ 

i=l i=l i=l + l 



Thus 



= & - Ax - E a.X.Zi - E Kp(K)P i z l - E Ap(A)f3iZi 
i=i i=i i=i+i 

k k I n 

= r - ^2a i X i z l - E Kp{K)PiZi - E X iP( X i)Pi z i~ E Pi A P( A )*i 

i=l i=l i=fc+l i=l+l 



and 



A a r m = A a r - ]T a i \" +1 z l - ]T X^p^i)^ 

i=l i=l 
I n 

- E K +1 p(*i)fozi - E 

z=fc+l z=i+l 
/c / n 

E ftA^, + e PiK* + E 

i=l i+1 i=/+l 



i=1 i=l 



- e \ ?+1 p(w*i - E -miM* 

j=fe+l i=J+l 



1=1 
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/ n 
i=k+l i=/+l 

Since DGMRES augmented with eigenvectors makes ||A a r m || minimized, we 
nave a, = — - — ^ a+i — ; — and 

l n 

A a r m = E (A A ? ~ A? +1 p(Ai)A)^ + E ^ " (7) 

From Lemma 1, (@|, and ([7]), we get the follow theorem. 

Theorem 3 Denote the spectrum of A by a(A) = {Ai, A2, • • • , Aj}. a a {A) :— 
{Afc+i, Afc+2, • • • ; A;} (fc < I). Choose tt a to be a closed domain in the X—plan 
that contains a a (A)\{0} but not A = 0, such that its boundary is twice differen- 
tiable with respect to arc-length. Denote by $ n (A) the conformal mapping of the 
exterior of £l a onto the exterior of the unit disk {w : \w\ > 1}. Then the vector 
x m generated by DGMRES augmented with eigenvectors satisfies 

\\A a r m \\<Km a+2 ^p™, (8) 

for all m, where K is a positive constant independent of to, k — max{kj : kj — 
ind(A- Xj),Xj e a a (A) \{0}} and p a = |l/$ o (0)| < 1. 

Compared © with ^ , we know the convergence of DGMRES augmented 
with eigenvectors is faster than that of DGMRES. 

The implementation of DGMRES augmented with eigenvectors is different 
from that of GMRES. We derive the algorithm step by step. 

1. Let the initial vector xq — 0. Compute (3 — ||^4 a ro|| and set v\ = 
{A a r )/f3. 

2. Orthogonalize the vectors A a ro, ■ ■ ■ , A m ~ 1 rQ via Arnoldi-Gram-Schmidt 
process. 

For j = 1,2, m — a. 
For i = 1,2, ---J 
Compute hij — (Avj,Vi). 
v~^i = Avj - htjVi 
end i 

Let h j+1 j = ||% || 
Set Vj+i ~ Vj/hj+x^ 
end j 

Consequently, we get a orthogonal matrix V m - a — [vi, V2, • ■ • , v m ^ a ] and a 
(m — a + 1) x (to — a) Hessenberg matrix H . 

3. Compute k approximated smallest magnitude nonzero eigenvalues of A 
and add their corresponding eigenvectors Z\ , • • • , z% to the subspace. Let = 
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[hij]{ m ~ a )x{m-a)- Obviously V^_ a AV m - a = HK When m < q - 1, where 
q is the degree of the minimal polynomial of A, the has full rank. We 
suppose Ai, A2, • • • , Afc are k smallest eigenvalues of and y%, j/2, • • • > Vk are 
the eigenvectors corresponding to them. Then := V m -aVi (i = 1,2, ••-,&) 
are the approximated eigenvectors of A. We add them to the subspace and 
orthogonalize. 

For j = m — a + 1 , • • • , m — a + fc 

Fori=l,2,---,j 

Compute ftij = ( J 4%-_ m+0 ,u i ). 

Compute %+i = Azj_ m+a - hijVi 

end i 

Set h j+ x,j = \\vj+i 
Set Uj+i = Vj+i/hj+ij 
end j. 

Denote = [vi,v 2 , ■ ■ ■ ,v m - a ,Zi, ■ ■ ■ ,z k ], = [hij] m - a +k+i,m-a+k, and 

7( ) = [vi,v 2 , ■ ■ ■ ,v m - a ,v m - a +i, ■ ■ -,v m - a +k+i]- It is easy to see AW = V^H W 
and x m = x + Wy m for some y m 6C" X ". It follows that r m = ro — AWy m = 

r Q -V^H i0) y m . So we get 

A a r m = A a r - A a+1 Wy m = pv x - A a+1 Wy m . (9) 
From this, we should continue to deal with A a+1 W / . 

A a+1 W = A a (AW) = A a V (0) H {0} 
= A a - 1 {AV^)H ( ° ) 

„a-lr lTrWTrW 
— A [Vi, V 2 , ■ ■ ■ , V m ^ a , U 7n _ a+1 , W m _ a+ 2, ' ' • , Wm-a+fc+2j-n n 

= A^V^H^H^ 



A a-(a-l)y(a-l)jj( a -l) . ..^(0) 



Denote H (a) ■ ■ ■ H {0) by H. From the above, we get 

A a+1 W = V (a) H. (10) 

Since r m is the residual vector produced by DGMRES augmented with eigen- 
vectros, from fl9l) and (1101) . we have 



\\A a r m \\ = (3 Vl - A a+1 Wy m \\ = - V^(H)y m \\ 
= ||/3ei - Hy m \\ = min||/3ei — Hy\\ 

Apply the QR factorization of H: if = QR, where R is upper triangular. Thus 
y m satisfies Ry m — /?(Q T ei) and x m follows. 
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4. The practical computation of H. From the above we know 



A[vi, ■ ■ ■ ,V m - a ,Zi, ■■■ ,Zk] = [Vi, ■ ■ ■ ,V m - a ,V m - a+ i, ■ ■ ■ ,V m - a+ k+l]H^°\ 



and 



A[vi,V2, • • • , V m - a , V-m-a+1, 1 - 1 , 

[vi,v-2,- 



,(!) 



,(!) 



where 



H 



(i) 



/ /ill 


hi2 ■ 


hl.m — a 


hl.m-a+l 


/ll,m-o+fc+l 


/i21 


h 2 2 ■ 




h2,m-a+l 


/ll,m-o+fc+l 





h 32 ■ 


^3,m-a 


/l3,m-o+l 


/l3,m-o+fc+l 





• 


fom—a-\-l,m—a 




^m-a+l,m-a+fc+l 





• 





hm—a+2,m—a+l 


^m-a+2,m-a+fc+l 


V o 


• 








• ^m-o+fe+2,m-a+fe+l / 



Tr(°) . T7(!) 



/lij is the new entries we need to compute when from H to H , that is, 

CC' 

Next 



we need to compute I 2 ( m a+2)+fc]x(fc+i) entr j eg additional. 



A[ Vl , V2, ■ ■ ■ , fm-oi "m-o+li u L-a+2i ' 



(1) i 



= [VI, 



lT 7(2) 

, fm-o+li w m-a+2i w m-a+3i ' ' ' ' W m-a+fc+2-l 



.,(!) 



„(2) 



„(2) 



where 



H 



(i) 





hi2 ■ 


hl.m-a+l 


hl.m-a+2 


hl.m-a+k+2 


h 2 i 


h 2 2 ■ 


h2,m-a+l 


h2,m-a+2 


tl2,m-a+k+2 





h 32 ■ 


^3,ro-a+l 


^3,m-a+2 


^3,m-a+fc+2 





• 


hm—a+2,m.—a+l 


h"m — a+2,m — a+2 


hm-a+2,m-a+k+2 


V o 


• 








■ h m - a +k+2,m-a+k+2 



Similarly, from H W to 7? (2) , there are [2(™-a)+5+fc]x(fc+i) need to be com _ 



(3) 

puted. Following this way, we can continue to get H , • 
computed. 

We summarize the above as the following algorithm. 



• , and H can be 



Algorithm 4 DGMRES augmented with eigenvectors 

1. Pick initial vector x n and compute r = b — Ax and (3 = ||A°ro||. 

2. Apply Arnoldi process to A a ro 7 • • • , A" l_1 ro. 
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vi = r /||r || 

For j = 1 , • • • , m — a 

Fori = l,---,j 

Compute hij — (Avj , Vi ) 

Compute Vj+i — Avj — hijVi 

end i 

Set h j+hj = \\vj\\ 
Vj+i =v j+1 /h j+ltj 
end j 

4- Denote [/iij]( m - a )x(m-a) by H and compute its k eigenvectors yi,---,Uk 
corresponding to k smallest magnitude eigenvalues. Compute Z{ = Vyi- We 
add z[s to the subspace and denote W — [v\, ■ ■ ■ , v m - a , z\, ■ ■ ■ , Zk]- We apply 

Arnoldi process to AW and get the matrix and . 
For j = m — a + 1 , • • ■ , m — a + k 
For i = 1,2 • • ■ ,j 
Compute h i: j = (A^_( m _ a ), Vi) 
Compute Vj+i = A^_( m _ a ) - h i3 Vi 
end i 

Set h j+hj = \\vj\\ 

Set v j+1 =v j+1 /h j+ltj 

end j 

5. Compute H {1) ,---,H {a) . 
Fort = l,2,---,a 

For j = m — a + t, ■ • ■ ,m — a + k + t 
Fori = 1,2, ---J 
Compute h^ — {Avj,Vi) 
Vj+i = Avj - h^Vi 
end i 

h j+i,i = \\ V A\ 
Vj+i =v j+1 /h j+1J 
end j 

Compute H = ^ • • • H ^ . 

6. Compute the QR factorization of H : H = QR, where R is upper triangular. 
Set c = l3Q T e\. Solve the least-square problem R m y m = fi{Q T e\) and x m = 
x + Wy m . If r m := b — Ax m satisfies the residual norm then stop; else set 
xq = x m and go to 2. 

4 Numerical examples 

For convenience, We denote the DGMRES augmented with eigenvectors by 
ADGMRES. In this section we present some numerical examples to compare 
ADGMRES to DGMRES. All the experiments were performed in MAT LAB® 
7.5 on an Inter Core 2 Duo 2000MHz PC with main memory 1000M. 
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Example 1 We take A to be a 12 by 12 singular matrix which has the 
following Jordan canonical form 
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The right side b = (1, 1, • • • , 1) T . The ADGMRES uses m = 6, k = 1 and DGM- 
RES uses m = 7 and denote them by DGMRES(6,1), DGMRES(7) respectively. 
So they use the same size subspace. When < e , wc stop. The conver- 

gence curve for ADGMRES and DGMRES are indicated in Fig 1. It is easy to 
see ADGMRES convergence faster than DGMRES 




iterations 



Figl. The convergence curves for ADGMRES(6,1) and DGMRES(7) 

Example 2. In this example the matrix A is the same as the above but 
let 07.7 = 1000. Obviously such A has a larger eigenvalue 1000. We expect 
ADGMRES(4,1) convergence much faster than DGMRES{h). As indicated in 
Fig 2., it is the case. 
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10" 

ADGMRES 

DGMRES 



10"' 



10 



1000 2000 3000 4000 5000 6000 7000 

iterations 

Fig2. The convergence curves for ADGMRES(4,1) and DGMRES(5) when A has a 

very lager eigenvalue. 

From Fig 2. we can see after 70 iterations the DGMRES(5) stagnates but 
after 5000 iterations the curve of ADGMRES(4,1) decreases steeply. 

Example 3. In this example we take 077 — 0.001, that is, A has a very 
small eigenvalue. From Fig 3. we see DGMRES(5) stagnates after 50 iterations 
but ADGMRES (4,1) still works well. 



10" 



10"' 



- ADGMRES 
- DGMRES 



10" 



o 10 



10" 



2000 4000 6000 8000 10000 12000 

iterations 



Fig3. The convergence curves for ADGMRES(4,1) and DGMRES(5) when A has a 

very small eigenvalue. 
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Example 4. This example comes from [5], 
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In [5] The author found that DGMRES(2) converges faster than DGM- 
RES(3). Since A has only one nonzero eigenvalue 1, we infer from Theorem that 
the convergence of ADGMRES(2,1) is as the same as that of DGMRES(2). From 
Table 1 we observe it is the case. At the same time we also see DGMRES(3) 
stagnates. The residual norm produced by the three methods is indicated in 
the following table. 



runs 


ADGMRES(2,1) 


DGMRES(3) 


DGMRES(2) 


50 


0.0038 


0.00279 


0.0038 


100 


1.23 x 10~ 5 


0.00276 


1.23 x 10~ 5 


200 


1.71 x 10~ 9 


0.00276 


1.72 x 10~ 9 


300 


6.155 x 10" 14 


0.00276 


6.145 x 10" 14 



Table 1. The residual norm produced by ADGMRES(2,1), DGMRES(3), 
DGMRES(2) after 50, 100, 200, 300 runs respectively. 

5 Conclusion. 

The DGMRES method augmented with eigenvectors can improve the conver- 
gence, especially when the matrix has small or large eigenvalues. The DGMRES 
method often stagnates (see [H[S]). The DGMRES augmented with eigenvec- 
tors is a good choice to overcome the stagnation. When k, the number of the 
eigenvectors added to the subspace, is large, the method is expensive. So in 
practical use we usually choose a small k. 
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